pacman::p_load(sf, raster, spatstat, tmap, tidyverse)Hands-on_Ex02
Installing and Loading the R packages
Spatial Data Wrangling
Importing the spatial data
childcare_sf <- st_read("/Users/sharon/OneDrive - Singapore Management University/ISSS626 Data/Hands-on_Ex/Hands-on_Ex02/geospatial/child-care-services-geojson.geojson") %>%
st_transform(crs = 3414)Reading layer `child-care-services-geojson' from data source
`/Users/sharon/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/ISSS626 Data/Hands-on_Ex/Hands-on_Ex02/geospatial/child-care-services-geojson.geojson'
using driver `GeoJSON'
Simple feature collection with 1925 features and 2 fields
Geometry type: POINT
Dimension: XYZ
Bounding box: xmin: 103.6878 ymin: 1.247759 xmax: 103.9897 ymax: 1.462134
z_range: zmin: 0 zmax: 0
Geodetic CRS: WGS 84sg_sf <- st_read(dsn = "/Users/sharon/OneDrive - Singapore Management University/ISSS626 Data/Hands-on_Ex/Hands-on_Ex02/geospatial/", layer="CostalOutline")Reading layer `CostalOutline' from data source
`/Users/sharon/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/ISSS626 Data/Hands-on_Ex/Hands-on_Ex02/geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 1 feature and 67 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 103.6091 ymin: 1.16639 xmax: 104.0858 ymax: 1.471388
Geodetic CRS: WGS 84mpsz_sf <- st_read(dsn = "/Users/sharon/OneDrive - Singapore Management University/ISSS626 Data/Hands-on_Ex/Hands-on_Ex02/geospatial/",
layer = "MP14_SUBZONE_WEB_PL")Reading layer `MP14_SUBZONE_WEB_PL' from data source
`/Users/sharon/Library/CloudStorage/OneDrive-SingaporeManagementUniversity/ISSS626 Data/Hands-on_Ex/Hands-on_Ex02/geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 323 features and 15 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
Projected CRS: SVY21Mapping the geospatial data sets
tmap_mode('view')tmap mode set to interactive viewingtm_shape(childcare_sf)+
tm_dots()tmap_mode('plot')tmap mode set to plottingGeospatial Data wrangling
childcare <- as_Spatial(childcare_sf)
mpsz <- as_Spatial(mpsz_sf)
sg <- as_Spatial(sg_sf)childcareclass : SpatialPointsDataFrame
features : 1925
extent : 11810.03, 45404.24, 25596.33, 49300.88 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
variables : 2
names : Name, Description
min values : kml_1, <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>100044</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>44, TELOK BLANGAH DRIVE, #01 - 19/51, SINGAPORE 100044</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td></td> </tr><tr bgcolor=""> <th>NAME</th> <td>PCF SPARKLETOTS PRESCHOOL @ TELOK BLANGAH BLK 44 (CC)</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>349C54F201805938</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20211201093837</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center>
max values : kml_999, <center><table><tr><th colspan='2' align='center'><em>Attributes</em></th></tr><tr bgcolor="#E3E3F3"> <th>ADDRESSBLOCKHOUSENUMBER</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSBUILDINGNAME</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSPOSTALCODE</th> <td>99982</td> </tr><tr bgcolor=""> <th>ADDRESSSTREETNAME</th> <td>35, ALLANBROOKE ROAD, SINGAPORE 099982</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSTYPE</th> <td></td> </tr><tr bgcolor=""> <th>DESCRIPTION</th> <td>Child Care Services</td> </tr><tr bgcolor="#E3E3F3"> <th>HYPERLINK</th> <td></td> </tr><tr bgcolor=""> <th>LANDXADDRESSPOINT</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>LANDYADDRESSPOINT</th> <td></td> </tr><tr bgcolor=""> <th>NAME</th> <td>ISLANDER PRE-SCHOOL PTE LTD</td> </tr><tr bgcolor="#E3E3F3"> <th>PHOTOURL</th> <td></td> </tr><tr bgcolor=""> <th>ADDRESSFLOORNUMBER</th> <td></td> </tr><tr bgcolor="#E3E3F3"> <th>INC_CRC</th> <td>4F63ACF93EFABE7F</td> </tr><tr bgcolor=""> <th>FMEL_UPD_D</th> <td>20211201093837</td> </tr><tr bgcolor="#E3E3F3"> <th>ADDRESSUNITNUMBER</th> <td></td> </tr></table></center> mpszclass : SpatialPolygonsDataFrame
features : 323
extent : 2667.538, 56396.44, 15748.72, 50256.33 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs
variables : 15
names : OBJECTID, SUBZONE_NO, SUBZONE_N, SUBZONE_C, CA_IND, PLN_AREA_N, PLN_AREA_C, REGION_N, REGION_C, INC_CRC, FMEL_UPD_D, X_ADDR, Y_ADDR, SHAPE_Leng, SHAPE_Area
min values : 1, 1, ADMIRALTY, AMSZ01, N, ANG MO KIO, AM, CENTRAL REGION, CR, 00F5E30B5C9B7AD8, 16409, 5092.8949, 19579.069, 871.554887798, 39437.9352703
max values : 323, 17, YUNNAN, YSSZ09, Y, YISHUN, YS, WEST REGION, WR, FFCCF172717C2EAF, 16409, 50424.7923, 49552.7904, 68083.9364708, 69748298.792 sgclass : SpatialPolygonsDataFrame
features : 1
extent : 103.6091, 104.0858, 1.16639, 1.471388 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
variables : 67
names : ID_0, ISO, NAME_ENGLI, NAME_ISO, NAME_FAO, NAME_LOCAL, NAME_OBSOL, NAME_VARIA, NAME_NONLA, NAME_FRENC, NAME_SPANI, NAME_RUSSI, NAME_ARABI, NAME_CHINE, WASPARTOF, ...
value : 205, SGP, Singapore, SINGAPORE, Singapore, Singapore, NA, NA, NA, Singapour, Singapur, СингапÑÑ, سÙغاÙÙرة, æ°å å¡, NA, ... Converting the Spatial* class into generic sp format
childcare_sp <- as(childcare, "SpatialPoints")
sg_sp <- as(sg, "SpatialPolygons")childcare_spclass : SpatialPoints
features : 1925
extent : 11810.03, 45404.24, 25596.33, 49300.88 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs sg_spclass : SpatialPolygons
features : 1
extent : 103.6091, 104.0858, 1.16639, 1.471388 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs Converting the generic sp format into spatstat’s ppp format
childcare_ppp <- as.ppp(childcare_sf)Warning in as.ppp.sf(childcare_sf): only first attribute column is used for
markschildcare_pppMarked planar point pattern: 1925 points
marks are of storage type 'character'
window: rectangle = [11810.03, 45404.24] x [25596.33, 49300.88] unitsplot(childcare_ppp)Warning in default.charmap(ntypes, chars): Too many types to display every type
as a different characterWarning: Only 10 out of 1925 symbols are shown in the symbol map
summary(childcare_ppp)Marked planar point pattern: 1925 points
Average intensity 2.417323e-06 points per square unit
Coordinates are given to 11 decimal places
marks are of type 'character'
Summary:
Length Class Mode
1925 character character
Window: rectangle = [11810.03, 45404.24] x [25596.33, 49300.88] units
(33590 x 23700 units)
Window area = 796335000 square unitsHandling duplicated points
any(duplicated(childcare_ppp))[1] FALSEmultiplicity(childcare_ppp) [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[38] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[75] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[112] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[149] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[186] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[223] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[260] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[297] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[334] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[371] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[408] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[445] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[482] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[519] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[556] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[593] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[630] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[667] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[704] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[741] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[778] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[815] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[852] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[889] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[926] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[963] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1000] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1037] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1074] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1111] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1148] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1185] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1222] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1259] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1296] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1333] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1370] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1407] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1444] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1481] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1518] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1555] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1592] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1629] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1666] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1703] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1740] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1777] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1814] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1851] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1888] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[1925] 1sum(multiplicity(childcare_ppp) > 1)[1] 0tmap_mode('view')tmap mode set to interactive viewingtm_shape(childcare) +
tm_dots(alpha=0.4,
size=0.05)tmap_mode('plot')tmap mode set to plottingchildcare_ppp_jit <- rjitter(childcare_ppp,
retry=TRUE,
nsim=1,
drop=TRUE)any(duplicated(childcare_ppp_jit))[1] FALSE# Transform the CRS to SVY21 (EPSG:3414)
sg_sf_projected <- st_transform(sg_sf, 3414)
# Check if the CRS transformation was successful
st_crs(sg_sf_projected)Coordinate Reference System:
User input: EPSG:3414
wkt:
PROJCRS["SVY21 / Singapore TM",
BASEGEOGCRS["SVY21",
DATUM["SVY21",
ELLIPSOID["WGS 84",6378137,298.257223563,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
ID["EPSG",4757]],
CONVERSION["Singapore Transverse Mercator",
METHOD["Transverse Mercator",
ID["EPSG",9807]],
PARAMETER["Latitude of natural origin",1.36666666666667,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8801]],
PARAMETER["Longitude of natural origin",103.833333333333,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8802]],
PARAMETER["Scale factor at natural origin",1,
SCALEUNIT["unity",1],
ID["EPSG",8805]],
PARAMETER["False easting",28001.642,
LENGTHUNIT["metre",1],
ID["EPSG",8806]],
PARAMETER["False northing",38744.572,
LENGTHUNIT["metre",1],
ID["EPSG",8807]]],
CS[Cartesian,2],
AXIS["northing (N)",north,
ORDER[1],
LENGTHUNIT["metre",1]],
AXIS["easting (E)",east,
ORDER[2],
LENGTHUNIT["metre",1]],
USAGE[
SCOPE["Cadastre, engineering survey, topographic mapping."],
AREA["Singapore - onshore and offshore."],
BBOX[1.13,103.59,1.47,104.07]],
ID["EPSG",3414]]Creating owin object
sg_owin <- as.owin(sg_sf_projected)plot(sg_owin)
summary(sg_owin)Window: polygonal boundary
34 separate polygons (no holes)
vertices area relative.area
polygon 1 4033 604202000.0 8.67e-01
polygon 2 38 414592.0 5.95e-04
polygon 3 16 19671.9 2.82e-05
polygon 4 537 24890800.0 3.57e-02
polygon 5 106 1949160.0 2.80e-03
polygon 6 261 10256000.0 1.47e-02
polygon 7 36 38729.5 5.56e-05
polygon 8 76 1087280.0 1.56e-03
polygon 9 62 920917.0 1.32e-03
polygon 10 96 573783.0 8.23e-04
polygon 11 34 126330.0 1.81e-04
polygon 12 13 20190.5 2.90e-05
polygon 13 18 34466.7 4.94e-05
polygon 14 37 239043.0 3.43e-04
polygon 15 320 4789260.0 6.87e-03
polygon 16 473 26400400.0 3.79e-02
polygon 17 21 47547.0 6.82e-05
polygon 18 77 1430070.0 2.05e-03
polygon 19 192 4796520.0 6.88e-03
polygon 20 132 2985250.0 4.28e-03
polygon 21 75 955310.0 1.37e-03
polygon 22 19 63805.9 9.15e-05
polygon 23 105 832107.0 1.19e-03
polygon 24 21 35903.4 5.15e-05
polygon 25 43 139236.0 2.00e-04
polygon 26 23 71722.4 1.03e-04
polygon 27 70 483000.0 6.93e-04
polygon 28 207 3885250.0 5.57e-03
polygon 29 16 35110.5 5.04e-05
polygon 30 26 72188.5 1.04e-04
polygon 31 156 2364690.0 3.39e-03
polygon 32 67 732165.0 1.05e-03
polygon 33 104 1100540.0 1.58e-03
polygon 34 82 1034550.0 1.48e-03
enclosing rectangle: [3040.59, 56097.76] x [16599.19, 50324.13] units
(53060 x 33720 units)
Window area = 697027000 square units
Fraction of frame area: 0.39Combining point events object and owin object
childcareSG_ppp = childcare_ppp[sg_owin]summary(childcareSG_ppp)Marked planar point pattern: 1924 points
Average intensity 2.760294e-06 points per square unit
Coordinates are given to 11 decimal places
marks are of type 'character'
Summary:
Length Class Mode
1924 character character
Window: polygonal boundary
34 separate polygons (no holes)
vertices area relative.area
polygon 1 4033 604202000.0 8.67e-01
polygon 2 38 414592.0 5.95e-04
polygon 3 16 19671.9 2.82e-05
polygon 4 537 24890800.0 3.57e-02
polygon 5 106 1949160.0 2.80e-03
polygon 6 261 10256000.0 1.47e-02
polygon 7 36 38729.5 5.56e-05
polygon 8 76 1087280.0 1.56e-03
polygon 9 62 920917.0 1.32e-03
polygon 10 96 573783.0 8.23e-04
polygon 11 34 126330.0 1.81e-04
polygon 12 13 20190.5 2.90e-05
polygon 13 18 34466.7 4.94e-05
polygon 14 37 239043.0 3.43e-04
polygon 15 320 4789260.0 6.87e-03
polygon 16 473 26400400.0 3.79e-02
polygon 17 21 47547.0 6.82e-05
polygon 18 77 1430070.0 2.05e-03
polygon 19 192 4796520.0 6.88e-03
polygon 20 132 2985250.0 4.28e-03
polygon 21 75 955310.0 1.37e-03
polygon 22 19 63805.9 9.15e-05
polygon 23 105 832107.0 1.19e-03
polygon 24 21 35903.4 5.15e-05
polygon 25 43 139236.0 2.00e-04
polygon 26 23 71722.4 1.03e-04
polygon 27 70 483000.0 6.93e-04
polygon 28 207 3885250.0 5.57e-03
polygon 29 16 35110.5 5.04e-05
polygon 30 26 72188.5 1.04e-04
polygon 31 156 2364690.0 3.39e-03
polygon 32 67 732165.0 1.05e-03
polygon 33 104 1100540.0 1.58e-03
polygon 34 82 1034550.0 1.48e-03
enclosing rectangle: [3040.59, 56097.76] x [16599.19, 50324.13] units
(53060 x 33720 units)
Window area = 697027000 square units
Fraction of frame area: 0.39First-order Spatial Point Patterns Analysis
kde_childcareSG_bw <- density(childcareSG_ppp,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian") plot(kde_childcareSG_bw)
bw <- bw.diggle(childcareSG_ppp)
bw sigma
305.2404 childcareSG_ppp.km <- rescale.ppp(childcareSG_ppp, 1000, "km")childcareSG_ppp.km <- rescale.ppp(childcareSG_ppp, 1000, "km")kde_childcareSG.bw <- density(childcareSG_ppp.km, sigma=bw.diggle, edge=TRUE, kernel="gaussian")
plot(kde_childcareSG.bw)
Working with different automatic badwidth methods
bw.CvL(childcareSG_ppp.km) sigma
4.53932 bw.scott(childcareSG_ppp.km) sigma.x sigma.y
2.160434 1.395302 bw.ppl(childcareSG_ppp.km) sigma
0.389485 bw.diggle(childcareSG_ppp.km) sigma
0.3052404 kde_childcareSG.ppl <- density(childcareSG_ppp.km,
sigma=bw.ppl,
edge=TRUE,
kernel="gaussian")
par(mfrow=c(1,2))
plot(kde_childcareSG.bw, main = "bw.diggle")
plot(kde_childcareSG.ppl, main = "bw.ppl")
Working with different kernel methods
par(mfrow=c(2,2))
plot(density(childcareSG_ppp.km,
sigma=bw.ppl,
edge=TRUE,
kernel="gaussian"),
main="Gaussian")
plot(density(childcareSG_ppp.km,
sigma=bw.ppl,
edge=TRUE,
kernel="epanechnikov"),
main="Epanechnikov")Warning in density.ppp(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE, :
Bandwidth selection will be based on Gaussian kernelplot(density(childcareSG_ppp.km,
sigma=bw.ppl,
edge=TRUE,
kernel="quartic"),
main="Quartic")Warning in density.ppp(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE, :
Bandwidth selection will be based on Gaussian kernelplot(density(childcareSG_ppp.km,
sigma=bw.ppl,
edge=TRUE,
kernel="disc"),
main="Disc")Warning in density.ppp(childcareSG_ppp.km, sigma = bw.ppl, edge = TRUE, :
Bandwidth selection will be based on Gaussian kernel
Fixed and Adaptive KDE
kde_childcareSG_600 <- density(childcareSG_ppp.km, sigma=0.6, edge=TRUE, kernel="gaussian")
plot(kde_childcareSG_600)
kde_childcareSG_adaptive <- adaptive.density(childcareSG_ppp.km, method="kernel")
plot(kde_childcareSG_adaptive)
par(mfrow=c(1,2))
plot(kde_childcareSG.bw, main = "Fixed bandwidth")
plot(kde_childcareSG_adaptive, main = "Adaptive bandwidth")
# Load the required libraries
library(spatstat.geom)
library(sp)
# Step 1: Extract pixel values (z) and grid coordinates (x, y)
x_coords <- kde_childcareSG.bw$xcol # x coordinates of the grid
y_coords <- kde_childcareSG.bw$yrow # y coordinates of the grid
z_values <- t(kde_childcareSG.bw$v) # Pixel values (transposed to match coordinate orientation)
# Step 2: Create a grid of the coordinates
grid <- expand.grid(x = x_coords, y = y_coords)
# Step 3: Create a SpatialPixelsDataFrame (used for gridded data)
sp_grid <- SpatialPixelsDataFrame(points = grid,
data = data.frame(value = as.vector(z_values)),
proj4string = CRS("+proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs"))
# Step 4: Convert to SpatialGridDataFrame
sp_grid_df <- as(sp_grid, "SpatialGridDataFrame")
# Check the structure
summary(sp_grid_df)Object of class SpatialGridDataFrame
Coordinates:
min max
x 3.040593 56.09776
y 16.599186 50.32413
Is projected: TRUE
proj4string :
[+proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1
+x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs]
Grid attributes:
cellcentre.offset cellsize cells.dim
x 3.247848 0.4145091 128
y 16.730924 0.2634761 128
Data attributes:
value
Min. : 0.000
1st Qu.: 0.000
Median : 0.109
Mean : 2.767
3rd Qu.: 4.598
Max. :35.129
NA's :9998 #gridded_kde_childcareSG_bw <- as.SpatialGridDataFrame.im(kde_childcareSG.bw)
#spplot(gridded_kde_childcareSG_bw)# Load the raster package
library(raster)kde_childcareSG_bw_raster <- raster(kde_childcareSG.bw)kde_childcareSG_bw_rasterclass : RasterLayer
dimensions : 128, 128, 16384 (nrow, ncol, ncell)
resolution : 0.4145091, 0.2634761 (x, y)
extent : 3.040593, 56.09776, 16.59919, 50.32413 (xmin, xmax, ymin, ymax)
crs : NA
source : memory
names : layer
values : -7.060279e-15, 35.12852 (min, max)projection(kde_childcareSG_bw_raster) <- CRS("+init=EPSG:3414")
kde_childcareSG_bw_rasterclass : RasterLayer
dimensions : 128, 128, 16384 (nrow, ncol, ncell)
resolution : 0.4145091, 0.2634761 (x, y)
extent : 3.040593, 56.09776, 16.59919, 50.32413 (xmin, xmax, ymin, ymax)
crs : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +units=m +no_defs
source : memory
names : layer
values : -7.060279e-15, 35.12852 (min, max)Visualising the output in tmap
# Load the tmap package
library(tmap)
tmap_mode("plot") # For interactive mapstmap mode set to plottingtm_shape(kde_childcareSG_bw_raster) +
tm_raster("layer", palette = "viridis") +
tm_layout(legend.position = c("right", "bottom"), frame = FALSE)
Comparing Spatial Point Patterns using KDE
# Load the dplyr package
library(dplyr)pg <- mpsz_sf %>%
filter(PLN_AREA_N == "PUNGGOL")
tm <- mpsz_sf %>%
filter(PLN_AREA_N == "TAMPINES")
ck <- mpsz_sf %>%
filter(PLN_AREA_N == "CHOA CHU KANG")
jw <- mpsz_sf %>%
filter(PLN_AREA_N == "JURONG WEST")pg <- mpsz_sf %>%
filter(PLN_AREA_N == "PUNGGOL")
tm <- mpsz_sf %>%
filter(PLN_AREA_N == "TAMPINES")
ck <- mpsz_sf %>%
filter(PLN_AREA_N == "CHOA CHU KANG")
jw <- mpsz_sf %>%
filter(PLN_AREA_N == "JURONG WEST")par(mfrow=c(2,2))
plot(pg, main = "Ponggol")Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all
plot(tm, main = "Tampines")Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all
plot(ck, main = "Choa Chu Kang")Warning: plotting the first 10 out of 15 attributes; use max.plot = 15 to plot
all
plot(jw, main = "Jurong West")Warning: plotting the first 9 out of 15 attributes; use max.plot = 15 to plot
all
# Load the spatstat.geom package
library(spatstat.geom)pg_owin = as.owin(pg)
tm_owin = as.owin(tm)
ck_owin = as.owin(ck)
jw_owin = as.owin(jw)childcare_pg_ppp = childcare_ppp_jit[pg_owin]
childcare_tm_ppp = childcare_ppp_jit[tm_owin]
childcare_ck_ppp = childcare_ppp_jit[ck_owin]
childcare_jw_ppp = childcare_ppp_jit[jw_owin]# Load the spatstat package
library(spatstat)childcare_pg_ppp.km <- rescale.ppp(childcare_pg_ppp, 1000, "km")
childcare_tm_ppp.km <- rescale.ppp(childcare_tm_ppp, 1000, "km")
childcare_ck_ppp.km <- rescale.ppp(childcare_ck_ppp, 1000, "km")
childcare_jw_ppp.km <- rescale.ppp(childcare_jw_ppp, 1000, "km")
# Calculate bandwidths using bw.ppl()
bw_pg <- bw.ppl(childcare_pg_ppp.km)
bw_tm <- bw.ppl(childcare_tm_ppp.km)
bw_ck <- bw.ppl(childcare_ck_ppp.km)
bw_jw <- bw.ppl(childcare_jw_ppp.km)
# Perform kernel density estimation with the calculated bandwidth
kde_pg <- density(childcare_pg_ppp.km, sigma = bw_pg, edge = TRUE, kernel = "gaussian")
kde_tm <- density(childcare_tm_ppp.km, sigma = bw_tm, edge = TRUE, kernel = "gaussian")
kde_ck <- density(childcare_ck_ppp.km, sigma = bw_ck, edge = TRUE, kernel = "gaussian")
kde_jw <- density(childcare_jw_ppp.km, sigma = bw_jw, edge = TRUE, kernel = "gaussian")
# Plot the kernel density estimates
par(mfrow = c(2, 2))
plot(kde_pg, main = "Punggol")
plot(kde_tm, main = "Tampines")
plot(kde_ck, main = "Choa Chu Kang")
plot(kde_jw, main = "Jurong West")
par(mfrow=c(2,2))
plot(density(childcare_pg_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian"),
main="Punggol")
plot(density(childcare_tm_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian"),
main="Tempines")
plot(density(childcare_ck_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian"),
main="Choa Chu Kang")
plot(density(childcare_jw_ppp.km,
sigma=bw.diggle,
edge=TRUE,
kernel="gaussian"),
main="JUrong West")
par(mfrow=c(2,2))
plot(density(childcare_ck_ppp.km,
sigma=0.25,
edge=TRUE,
kernel="gaussian"),
main="Chou Chu Kang")
plot(density(childcare_jw_ppp.km,
sigma=0.25,
edge=TRUE,
kernel="gaussian"),
main="JUrong West")
plot(density(childcare_pg_ppp.km,
sigma=0.25,
edge=TRUE,
kernel="gaussian"),
main="Punggol")
plot(density(childcare_tm_ppp.km,
sigma=0.25,
edge=TRUE,
kernel="gaussian"),
main="Tampines")
Nearest Neighbour Analysis
clarkevans.test(childcareSG_ppp,
correction="none",
clipregion="sg_owin",
alternative=c("clustered"),
nsim=99)
Clark-Evans test
No edge correction
Z-test
data: childcareSG_ppp
R = 0.52559, p-value < 2.2e-16
alternative hypothesis: clustered (R < 1)clarkevans.test(childcare_ck_ppp,
correction="none",
clipregion=NULL,
alternative=c("two.sided"),
nsim=999)
Clark-Evans test
No edge correction
Z-test
data: childcare_ck_ppp
R = 0.88703, p-value = 0.063
alternative hypothesis: two-sidedclarkevans.test(childcare_tm_ppp,
correction="none",
clipregion=NULL,
alternative=c("two.sided"),
nsim=999)
Clark-Evans test
No edge correction
Z-test
data: childcare_tm_ppp
R = 0.67231, p-value = 1.194e-11
alternative hypothesis: two-sided